An approximate solution to the Boltzmann equation for vibrated granular disks 
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The behaviour of the lower order moments of the velocity distribution function for a system of inelastic 
granular disks driven by vertical vibrations is studied using a kinetic theory. A perturbative kinetic theory for 
vibro-fluidised beds was proposed by Kumaran (JFM, v. 364, 163). A scheme to generalise this theory to higher 
orders in the moments is presented here. With such a method it is possible to obtain an analytical solution to the 
moments of the distribution function up to third order. 



I. INTRODUCTION 



II. BASIC FORMULATION 
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The dynamics of vibrated granular materials, its instabili- 
ties, and pattern formation are of some interest in the recent 
years as demonstrated by experiments of [[[]] and simulations 
of [^]. The theoretical description of such systems is compli- 
cated by the fact that it is a driven dissipative system charac- 
terised by highly inelastic collisions and hence the validity of 
equations of hydrodynamics is not clear at present [|[|. How- 
ever, it is possible to describe one idealised situation, where 
the dissipation due to particle collisions is small and the am- 
plitude of wall oscillations is small compared to the mean free 
path, as was shown in the kinetic theories [[|, ||]. Such a de- 
scription might be one of the starting points where we can as- 
certain with some confidence the rigour of the approach used. 
The present work is a continuation of such an approach. 

In this communication, we show that it is possible to ob- 
tain an analytical solution to the Boltzmann equation for a 
dilute bed of vibrated granular disks, by the method of mo- 
ments, correct upto third order in the moments of the distri- 
bution function. The method of approach followed holds the 
same principle as followed in [f|], except in the choice of the 
distribution function, which is done here by expanding it in 
the orthogonal set of Hermite polynomials. It is hoped that 
the analysis of this base state solution for its stability would 
give some clues to understanding the instabilities occuring in 
a vibro fluidised bed. 

In Sec. II, we present a general methodology to approx- 
imate the distribution function by expanding it in Hermite 
polynomials, and a procedure to solve the Boltzmann equa- 
tion. We present only the important results here, the reader 
is referred to [Q] for the details of the kinetic theory of vibro- 
fluidised beds. In Sec. |0l| we obtain an analytical solution 
when the formulation is restricted upto the third order in the 
moments. The results we obtain here are qualitatively similar 
to results of [Q], except that we have obtained an analytical 
solution whereas in the latter it was an approximate series so- 
lution. 
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The Boltzmann equation for the velocity distribution func- 
tion, f(x, u), for vertically vibrated beds is [|J: 



d t f + u*di*f - g 



where the collision integral is, 



df d c f 



dt 



(1) 



^ = a |du*dk(w*-k)(/{/^-/ 1 / 2 ) (2) 

To obtain an approximate solution to this equation we expand 
the distribution function about a Maxwell distribution in some 
space as: 



/(x,u) = ^/°[l + ^(x)^(u) 
Jo 



(3) 



where, /°(u) = j^e"* l T ° is the Maxwell distribution func- 
tion with the To left out of the definition for the simplifications 
to follow. The density p is also expanded about a leading order 
density field, 



P = Pa + 



(4) 



The functions ipj are chosen from a set of linearly indepen- 
dent function space. The parameters A/(x) are determined 
using the method of moments. In this method a set of func- 
tions, ipi(vL), equal in number to the number of unknowns are 
chosen and are multiplied with the Boltzmann equation, and 
integrated over the velocity space. This way we obtain a set 
of differential equations for the unknown parameters. 

The leading order density and temperature distribution were 
obtained in Ml for a dilute bed. 
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Here, N is the number of particles per unit width of the 
bed, g is the gravitational accelaration, e is the coefficient of 
restitution for particle-particle collisions, and (t/ 2 ) represents 
the mean square velocity of the vibrating surface. For sinu- 
soidal forcing with characteristic velocity Uq, this is given by 
(U 2 ) = U 2 /2. 
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The functions (pj are chosen from a set of linearly indepen- 
dent function space. The parameters Aj(x) are determined 
using the method of moments. In this method a set of func- 
tions, ipi(u), equal in number to the number of unknowns are 
chosen and are multiplied with the Boltzman equation, and in- 
tegrated over the velocity space. This way we obtain a set of 
differential equations for the unknown parameters. 

a. Non-dimensionalisation As a simplification we use 
the following non-dimensionalisation. u = u* / \/To, z — 
z*g/T . Substituting Eq. (||) in Eq. ([]]), multiplying by ipi 
and integrating over the velocity we obtain the steady state 
differential equation for the moments for variation only in the 
vertical direction, z as: 



ToVTo 



dui i> i u 7 .fd z p(l + Aj<pj) 



+ 7^ /dm f(i + A m 



dip* d c f 



du z 



dt 



(7) 



Here the second term has been simplified using the divergence 
theorem and the condition that the distribution function van- 
ishes for large velocities. When the collision term is integrated 
over the velocities, dui, it can be simplified from the form in 
Eq. (||) to an equivalent form (see in [Q], for example) 



/dniVi^ = y*duidu 2 dk (wk)p 2 

x + A mj )(l + A W2k ) Aipi (8) 

where, Aipi = [ip' u + ip' 2i — ipu — ip2i], is the total change 
in ipi due to collisions. The terms in Eq. (Q) and Eq. (^|) can 
be simplified by defining (.) = Jduf (.) and CM as the 
collision integral operator. Then we have from Eq. (Q), 



where, 
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CfA^i] with elastic collisions 

C[Aipi] inelastic collsions, excluding the above term 

C[AiPi tpu] + C[AiPi cp 2j ] 
C[Aipi <pij(p2k] 



Here, the superscript e for the collision terms above indicate 
that the collisions are considered to be elastic. 

A leading order equation is obtained by considering the 
Ajifj as perturbations to the Maxwell distribution and by con- 
sidering elastic collision in the collision integral. Thereby we 
obtain by setting Aj = 0, 
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(11) 



With, ipi — u z , the leading order density variation is given 
by d z p* = — 1, giving p* = —z + c or p = p°e~ z , where 
p° = Ng/T is the density at z = 0. Kumaran [f| had ob- 
tained the values of p° and To using a balance of the lead- 
ing order source and dissipation, for low densities are given in 
Eqs. (||) and (|(]). A high density correction to these values was 
obtained in [Q] in the leading order. In the present analysis we 
restrict ourselves to the low density limit. To obtain a first or- 
der balance in this limit, we neglect the quadratic term AjA k , 
and subtract out the leading order equation for low densities. 
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[{u z ipi) d z p + (u z ipiipj) (Ajd z p + pd z Aj)] 



gp 



ToV% 



du z 



, dipi 



= [C[AiPi] + Aj C[AiPi ipxj] + A k C[AiPi <p 2k ] 
v Jo 

+ AjA k C[Aipi<p 1 j<p 2 k}] (9) 

Further, dividing the above equation by g/To^/To and defin- 
ing p* = In p, we obtain the following equation for the un- 
known variables. 



Here we have used 1/e; = p°aTo/g = Na. Eqs. ( |12| ) is 
a set of coupled linear non-autonomous first order ordinary 
differential equations in the varibles Aj. If we incorporate 
the pertubation to density, Eq. (Q), in Eq. (^|), then the above 
equation for the first order quantities reads: 



Sijd z Aj + S?d zPl = 



Sij Gij 
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A. Boundary conditions 



cOq * , o I a a * i a a \ i nt> i n a The boundary conditions for the above equations may be 

obtained by the following method. A balance of the value of a 
P a To fQQe _|_ (jOi _|_ C^Aj + C 2e k AjA k ) (10) moment 4>i is considered when it collides with a wall moving 
9 with a velocity U. The change in the value of a moment due 
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to the collision is given by the relation: 
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Jdv* J d< /(u)^(u) = Jdu* z J d< /(u)^(u) 

— oo — oo 

U oo 

- |d<|d</(u)[0 4 (u')-^(u 



(14) 
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where the primed variable deotes the velocity of the particle 
after a collision. Simplifying the above equation we get, 



U oo 



oo oo 



du* z J d</(u)<Mu') = Jdu* z J d</(u)^(u) (15) 
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If the mean free path of the particles is large compared to the 
amplitude of vibration, then we can make an assumption that 
the wall is stationary at one position, then the above equation 
can be further simplified by averaging over different probable 
velocities of the bottom wall, which, in the case of a sine wave 
oscillation is, 



P(U)dU 



(16) 



if we write U = Uo sin#, and e = Uq/T , then the averaged 
equation in terms of the nondimensional quantities for the bal- 
ance of <pi, after making the substitution for the distribution 
function at z = 0, will be, 
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or simply as, 

E ij A j (0) = B i (18) 

Solving these simultaneous equations we obtain conditions 
satisfied by Aj at z = 0. 



III. THIRD-ORDER MOMENT FORMULATION 

We now take the specific case of the third order formula- 
tion. A third order formulation is the lowest order in which 
we can include the anisotropy in the distribution function. We 
now briefly explain the choice of the moments ipj. The dis- 
tribution function satisfies the following criteria, (i) the time 
averaged vertical flux is zero, (u z ) — 0, (ii) the distribution 
function is normalised to unity, (iii) it is symmetric in the hor- 
izontal velocities. It is convenient to choose the lower order 
moments from a multi-dimensional Hermite polynomials for 
the following reasons. These polynomials form a linearly in- 
dependent orthogonal basis, and the resulting equations are 
more convenient than a linearly independent set. In addition 
two of the above conditions will be automatically satisfied by 
the distribution function by setting the corresponding the cor- 
responding Aj to zero. Note that since the leading order dis- 
tribution function f° already satisfies the unit normalisation 



condition, we would require that Aj (ipj) = 0. The flux con- 



dition requires that Aj 



Z<Pj) 



0. Since 1 and u z are two of 



the linearly independent functions in the orthogonal Hermite 
polynomials, setting the coefficients of these two will ensure 
that the above two conditions will be satisfied at all orders of 
the polynomials. Thus the choice of the ipj becomes simple 
by incorporating the constraints on the distribution function 
directly. 

The set of multi-dimensional orthogonal polynomials can 
be obtained from the recurrence relation for the Hermite poly- 
nomials, 



H n+1 (x) = xH n (x) - nH n -\x) 
with H°(x) = l, H\x)=x 



(19) 



and the multidimensional set is then given by, 

H(u x ,u z ) = H m (u x )H n (u z ) even m, all n. (20) 

A symmetric distribution in the horizontal velocity can be en- 
sured by taking only even powers of u x , ie., even m in the 
above expression. The above polynomials are can be nor- 
malised to unity and the factor is 1 /ml n\. 

In the case of the third order approximation we obtain the 
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following polynomials, 

<pj = {-1 + u g , —3u z + u z ,—l + u x , —u z + u 2 u z }, 



(21) 



in which the first two members viz., 1 and u z have been omit- 
ted for reasons discussed above. We choose the same set for 
the moment generating functions ipi in Eq. (Oband the func- 
tions for the boundary conditions, fa in Eq. (|18[). This way 
we can get the same order of representation in the moment 
equations as well as in the boundary conditions. 

The moments of the distribution can be obtained by substi- 
tuting the expression (|T]) in Eq. (^J), which are: 



(ul) = (1 + 2A 3 ) 
(u 2 z ) = (l + 2A l ) 
(u 2 x u z ) = 2A 4 
(u 3 z ) = 6A 2 



IV. SOLUTION 
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With this simplification, the other variables can be written 
down in terms of A\ and its derivatives as 



M{z) 



-l 



24 



+ [-A! 1 (z) + 



A'{(z) 



A'{{z) A[ 4 \z) 



G 
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A A (z) = 



A 3 (z) = A 1 (z)-A>;(z)(^e : 
(-1 + e 2 ) 



-A[(z) 



A'l(z) 



A'({z) 



4%) 



(26a) 



(26b) 



(26c) 



With the moments considered in Eq. (|1]), the set of rela- 
tions in Eq. (12) for A^ can be written as 



d z A t 



LljAje- 



(23) 



where, 



To solve Eq. (|25|), we can make a reduction in order by the 
following substitution 



A"(z) = y(z) 



(27) 



rO q— I/O C< \ rl ^ Q— 1 /~ile il ^ q — 

^ij — D ik Wfcj ~ u kj), ^ij — —J % k ^tji °i — —^ik °fe ) 



With the transformation x = gT z , we obtain from Eq. (|25|), 
(Note: here x is just a transformation variable and not the 
cartesian co-ordinate) 



and by considering a moment ipi = u z in 
the density correction, 



Eq. (|3j) we have for 



d zPl + 2d z A 1 =2A 1 , (24) 
In the case of the third order approximation, we have 

L% = {{0, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}} 

4 = ir{{°> -I' °> §>> {-§» °> h °h {°. h °> -§}. 

{1,0,-1,0}} 
h \ = - Ii ^g^{0,l,0,3} 

We note here that the A, are independent of the density cor- 
rection pi in the first order approximation These equations can 
be rearranged into a single fourth order equation in A\, which 
is easily accomplished through a symbolic routine: 



A[ 4 \z) = -4A\*>(z)+[^ 



(3), 



-in 



2 -4j A'l(z) 
4(l-e 2 )e" 



(25) 



y - 



,y 



i 



(m 2 - A)y = Dx 2 



(28) 



where a dot accent denotes d x , and the constants are c = 
47r/e 2 , L> = -tt(1 - e 2 )/ef. Further with 



we get, 



w(x) — y{x)/x 2 , 



uH cw = D. 

x 



(29) 



(30) 



which is a modified Bessel equation of zeroth order. The so- 
lution to the homogeneous equation is given by: 



W h {x) = Ci I (-y/cx) + C 2 K (y/cx) 



(31) 



A particular solution to the inhomogeneous equation can be 
obtained by variation of parameters. The Wronskian of the 
above solutions can be written down from standard references, 



Ip(</cx) K {y/cx) 
lo(Vcx) K {^/cx) 



1 

x 



(32) 



The particular solution so obtained is given by: 

Dx 



w p (x) 



[io(Vcx) Ki(y/cx)+K (Vcx)h(^x)] . 

(33) 
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The most general solution to the inhomogeneous equation is 
then 



w(x) — Wh + w p (x). 



(34) 



The Eq. ( |27| ) can now be solved in terms of the independent 
variable x, by further reduction in order, by writing it, 



x (xAt(x) + Ai(x)) = y(x), 



B(x) = w(x) 



B(x) 



where, 



Ax(x) = B(x). 



(35) 



(36) 



(37) 



The solutions to these first order equations can be easily writ- 
ten as 



B(x) 



1 



dx' w(x') x + C3 



(38) 



and 



A x {x) = dx" 



I dx w(x ) x + C3 
The correction to the density is then from Eq. (Q), 

X 

Pi(x) = -2 Jdx' (Mp- + B(x'^+c 5 . 

The derivatives of A\ (z) can be written down as 

A' x {z) = -xA x (x) = -xB(x) 
A'l(z) = -x 2 w(x) 
A'l'(z) = -2x 2 w(x) - x 3 w(x) 



c 4 . (39) 



(40) 



(41a) 
(41b) 
(41c) 



These are used to evaluate the other unknown functions in 
Eqs. 



A. Boundary conditions 

Some of the constants of integration can be directly elimi- 
nated by requiring that the functions to take finite values for 
x = (or equivalently, z — > 00). Consider one of the linearly 
independent solutions of Wh(x) in E q. (pT| ), c 2 Kq(^/cx), 
which when integrated through Eqs. (38|) and (^) gives a 
term in A\ (x): c 2 K (\/cx)/c. Since Kq{^/cx) is singular as 
x — » 0, we require c 2 = 0. Similarly we require that C3 = 0, 
as this will lead to a singular term in A\ (x): C3 In x. 

The expressions in Eq. (33) cannot be integrated to obtain 
a closed form expression. Nevertheless, they can be easily 



evaluated numerically by converting the expressions Eq. 
and ( f39| ) to a definite integral, J^dx' added with some con- 
stant. The choice of the lower limit, xo, can, in general, be art- 
ibary to suit the convinience of matching the given boundary 
conditions; but such a choice, say J-,dx' , for the integral in 
Eq. (|8) will lead to singularities in the expression for A\(x) 
in Eq. (39|), similar to those obtained by the constant C3 term. 
To do away with this singularity, we choose xq = 0, expand 
the integral in Eq. ([38|) in Taylor's series, and subtract out the 
source of the singularity (which essentially is the constant of 
integration, i.e., the value of the integral evaluated at the lower 
limit). In expanding w(x) about x = 0, we note from Eq. ( p4| ) 
that w(0) — 0, therefore, 



B{x) 



dx' w(x') x 



dx (Wic 



w 2 x ' 





(W\ w 2 

— x H 

V 2 3 



w 2 2 



where w%, w%, . . . are constants coming from the Taylor's ex- 
pansion. Such an expansion is possible as the series converges 
for < x < 1. The above expression identically vanishes at 
x = 0, therefore we donot have to explicitly subtract out any 
singularity during numerical computation of the definite inte- 
gral. Such problems of subtracting out the singularity, how- 
ever, does not arise in the integral of Eq. ( [39] ) and any arbitary 
point, xq, can simply be chosen for the numerical integration, 
keeping in mind the boundary conditions. 

We are now left with two arbitary constants, c\ and C4. 
These are evaluated using boundary conditions in a manner 
similar to those used in [Q]. For the sake of convinience we 
choose the lower limit for the definite integral in Eq. ( [39| ) 
to be xq = 1, then C4 will simply b e eq ual to the value of 
A\(x = 1). It can be seen from Eq. ( |22b[ ) that A\ is directly 
propotional to the moment of the distribution function (u 2 ,), 
whose value at z = can be obtained by considering it as a 
vertical flux of momentum because of collisions with the wall. 
Assuming that the particles collide with the wall have a lead- 
ing distribution to be a Maxwellian, the flux of momentum 
along the vertical direction is given by, 

e 



<«2>U 



1 



substituting in Eq. (22b) we have, A\ 



l(x=l) 



C4 



(42) 
e/4, therefore, 
(43) 



Since, momentum is transferred only in the vertical direc- 
tion we have from Eqs. (E2I), A3| x= i = and At|a; = i = 0. 



The constant, ci, can now be evaluated from Eqs. (27), 
and (||). 



The constant A" (z)\ can be solved for from Eqs. 

4'(*)u=5 (^i+^)i 2=0 

= 4lf 



(44) 
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The other constant C5, is obtained by considering mass bal- 
ance for the density correction px- Since the total mass of 
the bed is balanced in the leading order density profile po, the 
balance for the correction to the density is given by 



dz* popi = 0, 



(45) 



or, 



dx p\ = 0. 



(46) 



We note here that we have not strictly incorporated the 
boundary conditions as discussed in Sec. II. This is because 
of a need to satisfy more strong boundary conditions of non- 
diverging solutions for z — > 00. The equations of the sort of 
Eq. dl8| ) would be useful in obtaining, for example, a series 
solution by numerical methods where we expand the solution 
in decaying functions such as Laguerre polynomials. 

Furthermore, the moments related to the horizontal direc- 
tion are not exactly satisfied, i.e., the boundary conditions re- 
lated to A3 and A4. While the value of A3 is not satisfied 
independently, but only partly in Eq. (|4|), A4 is never used. 
These can be used only in a higher order polynomial approxi- 
mation. Strictly therefore, only A\ is satisfied exactly. 
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FIG. 1: Density profile for Na = 3, e = 0.3. The analytical solu- 
tion obtained in this paper and the series solution of [[ij] show good 
agreement. There is a negative correction to the density at the bottom 
after which it decays exponentially as the leading order values. 



V. DISCUSSION 

The results obtained here are qualitatively the same as ob- 
tained in [Q]. At the time of the previous work, the only re- 
sults available were some experimental measurements of [^[j 
and only a qualitative comparison was possible to ascertain 
the validity of the theory. We have now compared the results 
of with a numerical Event Driven (ED) simulation of vibrated 
hard disks. The ED simulation was done with periodic bound- 
aries and with an approximate representation of the bottom 
wall The following figures show that the theory is in- 
deed good in the limit of its validity. The figures also show 
a comparison with the approximate series solution which was 
obtained in [Q]. In most cases the series solution provides a 
good approximate. 

The parameters chosen for the simulation are in correspon- 
dence with the limit of validity of the theory: e <C 1, and 
Na ~ O (1), (see ||], f° r a discussion). There is a small 
negative correction to the density near the bottom wall, as 
shown in Fig. |l|, due to the energy flow near the bottom wall, 
after which it falls of exponentially as the leading order den- 
sity. 

One important difference, in the formulation of the den- 
sity, between the present and the previous work is the fol- 
lowing. In the previous work a correction to the distribution 
function due to variations in the distribution function over dis- 
tances comparable to the particle diameter were described by 
using a small parameter, to- It was shown in a later work 
that this correction is essentially equivalent to the high density 



correction obtained on the lines of Enskog correction to dense 
gases, therefore it has been omitted in the present considera- 
tion of dilute bed. 

Due to the anisotropic nature of the energy input to the 
vibro-fluidised bed, the lower order moments of the distribu- 
tion clearly show the anisotropy. Figs. || and|| show the hor- 
izontal and vertical temperature profiles, while Fig. Q shows 
the vertical flow of energy. Anisotropies were also observed 
by us in deep bed simulations of disks which had the wave- 
like surface patterns although the nature of anisotropy was 
more pronounced even in the shape of the distribution func- 
tion itself (the vertical distribution had double peaks and the 
horizontal distribution had single peak and exponential tails). 
Could the presence of anisotropy be an important feature giv- 
ing rise to an instability in one direction ? A stability analysis 
of the solution from the present analysis model might help 
resolve this. The usual models based on hydrodynamic equa- 
tions do not take into account this anisotropy. 

As pointed out earlier, the density correction considered 
here is different from the one used in the previous work, in 
that the high density correction is not considered in first order 
in the distribution function. This results in a better prediction 
of the vertical temperature Fig. || and flux of energy Fig. ^ 
even when the density prediction is not expected to be good. 

To conclude, we have shown that (a) It is possible to obtain 
an analytical solution to the Boltzmann equation for vibro flu- 
idised bed in the low density limit correct upto the third order 
in the moments of the distribution function, (b) The quali- 
tative nature of results obtained here are similar to those ob- 
tained in [Q]. (c) The correction to the distribution function 
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200 




due to spatial variation of the order of a particle diameter was 
neglected here as this is turns out to be a a correction to a 
higher order in density ^ . With this it is seen that the den- 
sity still shows a negative correction at the bottom wall due to 
the energy flow, (d) The boundary conditions are overspeci- 
fied in the problem and we chose to satisfy exactly, only the 
ones involved in the momentum transfer the vertical direction, 
(e) Even with this restricted choice, the theoretical values for 
the different moments compare reasonably well with the sim- 
ulation particularly in the anisotropy exhibited, (f) This gives 
us some confidence to explore further with the stability of the 
solution, and with an higher order approximation to the distri- 
bution function if required. 
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FIG. 2: Horizontal temperature profile for No = 3, e = 0.3. A 
magnitude comparison with Fig. pi shows the degree of anisotropy in 
a vibrated bed. 
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FIG. 3: Vertical temperature profile for Na = 3, e = 0.3. The theo- 
retical values are closer to the simulation here, than for the horizontal 
temperature because of the boundary conditions imposed. 
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FIG. 4: Flux of energy in the vertical direction for Na = 3, e = 0.3. 
Here again, the theoretical values compare well in order of magni- 
tude, because of boundary conditions are imposed in this direction. 
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FIG. 6: Flux of energy prediction is good even when the density is 
high (z/~ 0.3). 



